require('igraph')
## Loading required package: igraph
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
require('ggplot2')
## Loading required package: ggplot2
require('reshape')
## Loading required package: reshape
require('gridExtra')
## Loading required package: gridExtra
setwd("~/git/ABadea/")

listGs<- list.files(path = "./graphml/", pattern = "*.graphml")

#read in covariates and graph list
#find those with common ids, sort by id

covariates<- read.csv("./covariates/predictors.csv",stringsAsFactors = F)
ids <- unlist( lapply(listGs,function(x)strtrim(x,6)))
common_id<- intersect(covariates$RUNNO , ids)

covariates <- covariates[covariates$RUNNO%in%common_id,]
covariates <- covariates[order(covariates$RUNNO),]  

listGs<- listGs[ids%in%common_id]
listGs<- listGs[order(listGs)]

#load and sort the tensor decomp data
load("./tensorDecompBinary.Rda")
load("./listData.Rda")
ids2 <- unlist( lapply(dataList,function(x)strtrim(x$id,6)))
diagC<- tensorDecomp$C[ids2%in% common_id,]
diagC<- diagC[order(ids2[ids2%in% common_id]),]



graphList<- lapply(listGs, function(x){
  read.graph( file = paste("./graphml/",x,sep = ""),format = "graphml")
})

AdjacencyList<- lapply(graphList, function(x){
  get.adjacency(x)
})

HemisphereList<- lapply(graphList, function(x){
  get.vertex.attribute(x,name="hemisphere")
})
colnames(AdjacencyList[[1]])
## NULL
order_by_hemi<- order(HemisphereList[[1]])


AbyHemiSphere<- (as.matrix(AdjacencyList[[1]]))[order_by_hemi,order_by_hemi]

plotHeatmap<-function(denseA, zlimMax=1){
  rownames(denseA) <- c(1:nrow(denseA))
  colnames(denseA) <- c(1:nrow(denseA))
  m<-melt(denseA)
  n<- nrow(denseA)
  p <- ggplot(m, aes(X1, X2)) +
    geom_tile(aes(fill =  value), colour = "white") + 
    scale_fill_gradient(low = "white",  high = "red",limits=c(0,zlimMax)) +
    theme(axis.ticks = element_blank(), axis.text.x = element_blank(),axis.text.y = element_blank()) 
  
  p
}

Adjacency Matrix

One random subject:

The adjacency matrix of the whole brain.

Segments: Bottom left: Left Hemisphere Upper right: Right Hemisphere

wholeBrain<- plotHeatmap(AbyHemiSphere)
wholeBrain+ geom_hline(yintercept = nrow(AbyHemiSphere)/2 ) + geom_vline(xintercept = nrow(AbyHemiSphere)/2 )

On average:

The adjacency matrix of the whole brain.

Segments: Bottom left: Left Hemisphere Upper right: Right Hemisphere

m<- length(AdjacencyList)


n<- nrow(AdjacencyList[[1]])

computeAvgA<- function(l){
  sumAdj<- matrix(0,n,n)
  for(i in 1:length(l)){
  sumAdj<- sumAdj + (AdjacencyList[[i]])[order_by_hemi,order_by_hemi]
  }
  as.matrix(sumAdj/m)
}

averageAdjacency<- computeAvgA(AdjacencyList)


avgWholeBrain<- plotHeatmap(averageAdjacency)
avgWholeBrain+ geom_hline(yintercept = nrow(AbyHemiSphere)/2 ) + geom_vline(xintercept = nrow(AbyHemiSphere)/2 )

Now we focus on the average graph:

lDegree<- rowSums(averageAdjacency[1:(n/2),1:(n/2)])
rDegree<- rowSums(averageAdjacency[(n/2+1):(n),(n/2+1):(n)])

df<- data.frame("ROI_Index"= rep(c(1:(n/2)),2),
                "Hemisphere"=as.factor(rep(c("L","R"),each=n/2)),
                "Degree" = c(lDegree,rDegree)
)

Degree plot by ROI index

ggplot(data=df,
       aes(x=ROI_Index, y=Degree, colour=Hemisphere)) +
  geom_line()

Degree distribution

ggplot(df, aes(x=Degree, fill=Hemisphere)) + geom_density(alpha = 0.7, bw = 3)

Degree distribution of the 3 age groups:

Young mice have higher level of connectivy (degree) than the middle and old.

df<- data.frame("ROI_Index"= rep(c(1:(n/2)),3),
                "Age"=as.factor(rep(c("Young","Middle","Old"),each=n/2)),
                "Degree" = c(rowSums(youngAvgA),rowSums(middleageAvgA),rowSums(oldAvgA))
)

ggplot(df, aes(x=Degree, fill=Age)) + geom_density(alpha = 0.7, bw = 3)

After dimension reduction:

Pairs plot of the first 5 dimensions:

require('GGally')
## Loading required package: GGally
ageGroup<-  (weeks<=63)*1 + (weeks>63 & weeks<=78)*2+(weeks>78)*3
ageGroup[ageGroup==1] = "Young"
ageGroup[ageGroup==2] = "Middle"
ageGroup[ageGroup==3] = "Old"

df<- data.frame("Latent Dimension Idx"= rep(c(1:30),m),
                "Age"=as.factor(rep(ageGroup,each=30)),
                "Latent Coordinate" = c(t(diagC))
)



K<- 5

seleDf<- diagC[,1:K]


# seleDf<- cbind((seleDf))

colnames(seleDf)<- c(  sapply(c(1:K),function(x){paste("x",x,sep = ".")}))

newDf<- as.data.frame(seleDf)
newDf<- cbind(newDf, ageGroup)
ggpairs(newDf,aes(col=ageGroup,alpha=0.4, binwidth=1))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Average scree plots of the 3 groups:

C1<- diagC[ageGroup=="Young",]
C2<- diagC[ageGroup=="Middle",]
C3<- diagC[ageGroup=="Old",]

df<- data.frame("Latent.Dim"= rep(c(1:(30)),3),
                "AgeGroup"=as.factor(rep(c("Young","Middle","Old"),each=30)),
                "Latent.Coor" = c(colMeans(C1),colMeans(C2),colMeans(C3))
)

ggplot(data=df,
       aes(x=Latent.Dim, y=Latent.Coor, colour=AgeGroup)) +
  geom_line()

From the low dimension representation, old group looks obviously distinct from the middle and young groups.